Spinodal phase separation in relativistic nuclear collisions 
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The spinodal amplification of density fluctuations is treated perturbatively within dissipative fluid 
dynamics for the purpose of elucidating the prospects for this mechanism to cause a phase separation 
to occur during a relativistic nuclear collision. The present study includes not only viscosity but 
also heat conduction (whose effect on the growth rates is of comparable magnitude but opposite), as 
well as a gradient term in the local pressure, and the corresponding dispersion relation for collective 
modes in bulk matter is derived from relativistic fluid dynamics. A suitable two-phase equation 
of state is obtained by interpolation between a hadronic gas and a quark-gluon plasma, while the 
transport coefficients are approximated by simple parametrizations that are suitable at any degree 
of net baryon density. We calculate the degree of spinodal amplification occurring along specific 
dynamical phase trajectories characteristic of nuclear collision at various energies. The results bring 
out the important fact that the prospects for spinodal phase separation to occur can be greatly 
enhanced by careful tuning of the collision energy to ensure that the thermodynamic conditions 
associated with the maximum compression lie inside the region of spinodal instability. 

PACS numbers: 25.75.-q, 81.30.Dz, 64.75. Gh, 64.60.an 
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I. INTRODUCTION 

It is expected that the confined and deconfined phases 
of strongly interacting matter may coexist at net baryon 
densities above a certain critical value and significant ex- 
perimental efforts are underway to search for evidence of 
the associated first-order phase transition and its critical 
end point: a systematic beam-energy scan is currently 
being performed at RHIC (BNL) to look for the critical 
point [1]; the CBM experiment at FAIR (GSI) will study 
baryon-dense matter and search for the phase transition 
; and the proposed NICA ( JINR) aims at exploring the 
mixed phase Q. 

These studies are rather challenging, not only because 
it is inherently difficult to extract the properties of equi- 
librium bulk matter from collision experiments, but also 
because there is currently no suitable transport model 
available for guiding these efforts. As a result, there is 
currently an urgent need for identifying experimentally 
observable signals of the phase structure. 

The present paper focusses on the possibility that the 
mechanism of spinodal phase decomposition may have ef- 
fects that could be exploited as signals of the phase tran- 
sition. Spinodal decomposition is a well-known generic 
phenomenon associated with first-order phase transitions 
that has been studied in a variety of substances and also 
found industrial application [j|. Furthermore, nuclear 
spinodal fragmentation [5| was observed in nuclear col- 
lisions at intermediate energies [6| several years ago. A 
preliminary study 7] found grounds for guarded opti- 
mism that spinodal separation between the confined and 
deconfined phases could in fact occur in relativistic colli- 
sions and we have therefore undertaken the present more 
refined analysis. 



schematic equation of state based on a generalized clas- 
sical gas, the present uses a more realistic equation of 
state obtained by interpolating between a hadron gas and 
a quark-gluon plasma. An advantage of this procedure 
is that it automatically incorporates the increase in the 
number of degrees of freedom in the dense (deconfined) 
phase, a peculiar but important characteristic of strongly 
interacting matter. Building on the developments in Ref. 
[T] , we take account of finite-range effects by including a 
gradient term in the equation of state. This refinement 
is essential for obtaining a physically meaningful descrip- 
tion since it ensures both that there is an interface tension 
between the two coexisting phases and that the spinodal 
growth rates subside at large wave numbers. 

We again carry out our studies within the framework 
of fluid dynamics, because this type of transport descrip- 
tion has the distinct advantage that the complicated and 
still poorly understood microstructure of the system en- 
ters only via the equation of state and the transport co- 
efficients. A general discussion of the fiuid-dynamical 
description of first-order phase transitions was given re- 
cently in Ref. [1]. 

Although the dispersion equation in [7| was derived 
with both shear and bulk viscosity included, the actual 
calculations were done for ideal fluid dynamics. In the 
present study, the dynamical calculations include not 
only viscosity but also heat conductivity, which proves 
to be as important as viscosity while affecting the spin- 
odal growth rates oppositely, as we demonstrate quan- 
titatively. The associated cubic dispersion equation is 
derived directly from relativistic fluid dynamics. The 
medium dependence of the transport coefficients is ex- 
pressed in terms of the local values of the enthalpy den- 
sity and the particle spacing, an approximation that ap- 
plies not only in the baryon-poor regime but also in the 



The strength of the transport coefRcients for arbitrary 
density and temperature can thus be related to the val- 
ues obtained from analyses of the RHIC data. 

We seek to construct plausible dynamical phase trajec- 
tories by invoking results from earlier calculations with 
various transport models and we examine in particular 
the crucial importance of using a collision energy for 
which the maximum compression occurs inside the spin- 
odal phase region. Once the phase trajectory of the colli- 
sion system has been specified, we may integrate the spin- 
odal growth rate along the dynamical history and thus 
calculate the resulting degree of amplification as a func- 
tion of the wave number of the density perturbation. We 
do that for a range of dissipation strengths that bracket 
those expected from the analysis of the RHIC data. 

The present, more refined, studies suggest that spin- 
odal phase decomposition may indeed be triggered during 
nuclear collisions within a certain (likely relatively nar- 
row) optimal energy range. This expectation is rather 
insensitive to the still poorly known strength of the trans- 
port coefficients. Such a spinodal phase separation would 
result in an assembly of plasma drops embedded in a 
hadron gas and our present analysis permits us to esti- 
mate the typical drop size. 

The presentation is organized as follows. We first 
discuss the expected thermodynamic phase structure of 
strongly interacting matter within the framework of the 
specific equation of state that we have constructed. We 
then turn to dissipative fiuid dynamics within which we 
derive the dispersion equation for the collective modes in 
bulk matter. Subsequently we develop expressions for the 
transport coefficients in baryon-rich matter and use those 
in calculations of the spinodal growth rates. Finally, we 
obtain quantitative results for the degree of spinodal am- 
plification experienced by the bulk of the collision system 
as it evolves along various plausible phase trajectories. 
The construction of the equation of state and the associ- 
ated spline procedure are described in appendices. 
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FIG. 1: The {p,T) phase diagram for bulk matter showing 
the phase coexistence boundary (outer solid), the isothermal 
spinodal where vt = (dashed), and the isentropic spinodal 
where Hs = (inner solid), together with the critical point 
(for details, see App. |X]). The dispersion relation shown in 
Fig. [2] was calculated at the square, while the results shown 
in Fig. [3]were obtained along the two straight lines. 



of the system is characterized by the temperature T and 
the (net) baryon density p = ub — ng. (In the phase re- 
gion of primary interest, the temperature is relatively low 
and the chemical potential relatively high, so the equilib- 
rium population of antibaryons is relatively insignificant, 
rig <^ riB-) The key thermodynamic function is then the 
free energy density,s frip), from which the other quanti- 
ties may be obtained. 

Chemical potential : urip) — dpfrip) , (1) 

Pressure : pt{p) = pdpfrip) ~ hip) , (2) 

Entropy density : (Tt{p) = -drfTip) , (3) 

Energy density : erip) = frip) ~ TdTfrip) ■ {^) 

We may also express the isothermal sound speed vt, 



II. THERMODYNAMIC PHASE STRUCTURE 



dp 



P'n2, 



^T-^Tzifp) - irJ^p-^p) = t-J'Mp) ■ (5) 



In order to make quantitative studies, we need to em- 
ploy a specific equation of state that is plausibly realistic 
and, in particular, exhibits the expected phase structure. 

Although significant progress has been made in un- 
derstanding the thermodynamic properties of the con- 
fined and deconfined phases separately, our current un- 
derstanding of the phase coexistence region is not yet on 
firm ground. We therefore employ a conceptually sim- 
ple approximate equation of state in which the region of 
phase transformation is described by means of a suitable 
interpolation between an idealized hadron gas and an ide- 
alized quark-gluon plasma. The details of this construc- 
tion are described in Appendix [^ while the resulting 
phase structure is shown in Fig. [TJ 

For the present study, it is convenient to work in the 

rj^.Tinnifj^l renrpsentj^.tinn "u/hpre thp therTnndvna.mir statp 



where h^^ (p) — px (p) + £t {p) is the enthalpy density, as 
well as the isentropic sound speed Vs, 
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where the specific heat is c-u = drSr^p) = TdrcrTip)- 

Different manifestations of the system that have the 
same values of temperature, chemical potential, and 
pressure are in mutual thermodynamic equilibrium and 
may thus coexist. As Fig. [1] shows, such phase coexis- 
tence occurs for temperatures below the critical value, 
T <Tc~ 142.5 MeV. The corresponding density values, 
p^ and pf", are traced out (outer solid curve). Inside 
this boundary, it is thermodynamically favorable for uni- 
form matter tn spnaratp into thp tw^n nnPYistincr nhasps 



However, in the part of the phase coexistence region that 
is close to the boundary any small deviation from uni- 
formity is thermodynamically unfavorable and here the 
system is mechanically metastable, as signalled by the 
fact that w|n > in this phase region. 

This situation changes radically at the spinodal bound- 
ary where the isothermal sound speed vanishes, vt = 
(the dashed curve in Fig. [l]): Inside this boundary even 
infinitesimal deviations from uniformity are thermody- 
namcially favorable and bulk matter is thus mechanically 
unstable. As a consequence, small density fluctuations 
may be amplified and thereby cause the system to un- 
dergo a spontaneous phase separation. The phase region 
of spinodal instability extends in temperature all the way 
from zero up to Tc- 

Finally, Fig. [T] also shows the boundary where the isen- 
tropic sound speed vanishes, Vg = 0. Inside this smaller 
phase region, the spinodal amplification process can oc- 
cur without any entropy production. This special region 
does not extend all the way up to Tc and it narrows con- 
siderably as the temperature approaches its upper bound. 

The equation of state applies to idealized uniform mat- 
ter. However, since we are interested in following the 
evolution of density disturbances, it is essential to in- 
clude finite-range effects. Indeed, the physical coexis- 
tence between two different phases along a common inter- 
face could not be realized without taking proper account 
of the density gradients, nor could the associated inter- 
face tension be obtained. We shall therefore augment the 
bulk thermodynamics with a gradient term as proposed 
in Ref. [3| (see Sect. lIIIC 1|) . which then extends the va- 
lidity of the equation of state to non-uniform systems. In 
particular, it is possible to describe the diffuse interface 
between two coexisting phases and the associated ten- 
sion. Furthermore, as we shall see, the gradient term is 
essential for obtaining a physically reasonable dispersion 
relation because it stabilizes disturbances having a small 
spatial scale. 



III. DISSIPATIVE FLUID DYNAMICS 

We wish to employ dissipative fluid dynamics for our 
dynamical studies. Fluid dynamics is convenient because 
the specific miscroscopic structure of the matter under 
consideration enters only via the equation of state and a 
few transport coefficients. On the other hand, the treat- 
ment relies on the assumption of approximate local equi- 
librium which may generally be questionable in nuclear 
collisions. Fortunately, our applications are to collisions 
of relatively modest energies and, moreover, to the later 
stages in the evolution. Thus the conditions for applica- 
bility should be reasonably favorable. 

We base our treatment on the relativistic formulation 
by Muronga [y]. The four- velocity of the local flow is 
'"'' — ill 1'^) Ei-nd the symmetric tensor A'^'^ = g^'^ — u^^u" 
projects onto the 3-space orthogonal to -u^, A'^'^Ui, = 0. 

TVip snarp-timp dprivativp Hprmnnnsps. r)'' = iit^D -\- V^. 



where the convective time derivative is D = u^^d^ = 
7[9t -I- V • V] while the gradient is V^ = A^^c^i,. 

The equations of motion simplify when expressed in 
a specific reference frame. For the present study, it is 
convenient to use the Eckart frame, which is defined in 
terms of the charge flow and is usually employed in non- 
relativistic scenarios. Then, since the local charge flow, 
y^, vanishes by definition, the charge four-current den- 
sity is N^^ — pu^, where p is the charge density in the 
local flow frame. However, it may generally be preferable 
to use the Landau frame because it is defined in terms of 
the energy flow and is thus meaningful also for chargeless 
fluids and fluids with several conserved charges for which 
there is no unique generalization of the Eckart frame. 



A. Small disturbances 

We consider the early evolution of small deviations 
from uniformity. We generally assume that these are 
planar and harmonic, p(r,t) = po + pkexp{ikx — iujt) 
with 5 = \pk\/ Po <SC 1, and similarly for the other quan- 
tities. We may then ignore terms of order 0{5'^) and 
higher. It follows that the associated flow velocities are 
small, w <C 1 since 0{v) = 0{8), and thus we have 



I'-- 



-1/2 



l + 0{5^)Kl. 



Generally the energy flow is given by W^ 



-hV 



q^ , where q^ is the heat fiow and h — p + e the enthalpy 
density. Since the local charge fiow V^ vanishes when we 
use the Eckart frame (see above), the energy fiow equals 
the heat flow, W^" = g^, and is 0{5). Hence 0{Wu) = 
0{5'^), and the energy- momentum tensor simplifies to 



T^"" = £u^u''-pA^ 



HA^ 



(7) 



Here e = u^T^^'^Uy is the energy density in the local fiow 
frame and p -I- H = — ■iA^^T'^'^ is the sum of the local 
isotropic pressure p and the pressure induced by the bulk 
viscosity which enters through the bulk pressure. 



H 



<yuu^ 



<y^v' 



^coy 



-CV-v. (8) 



Furthermore, the heat flow is q^ = u^T^'^'A^, while the 
shear viscosity enters via the stress tensor 

TT^" = r;[A^A^ + A^A^ - fA^'^A^.lV'u" (9) 
« r;[Af A^^ + A^AJ- - fA'^-AylV^z;^' , (10) 

where we have used that only the spatial components 
of V^ contribute to leading order in d and, moreover, 
any derivatives of it^ w 1 can be ignored so only the 
spatial components w* ~ v' contribute. Furthermore, 
since A° = w*, which is 0{6), only the spatial elements 
A* « 6ij contribute. Hence only the 3x3 spatial part tt 
is non- vanishing. It has the following elements, 

n'^ « -4d,v' + dy - U.^dkv"] . (11) 



Thus, for small deviations from uniformity, the spatial 
part T of the energy-momentum tensor is given by 



If the spatial variation of the viscosity coefEcients 77 
(shear) and C (bulk) may be ignored, we have 



V T 



Vp - r]Av ~ [\ri + C] V(V • v) , 



(13) 



The gradient simplifies further in a semi-infinite geome- 
try (where the only spatial variation is in the x direction), 

d,p-[h + (:]div. (14) 



V T 



B T 



Thus only the effective viscosity ^ = |?7 + C enters. It 
follows that a uniform stretching, v{x) ~ x, is dissipation 
free. We also wish to point out that the shear viscosity 
contributes even though the flow has no shear. 

It is interesting to note that the above result reflects a 
general feature of isotropic expansions in N dimensions. 
To see this, assume that p{r) — p{r) and v{r) — v{r)f. 
The viscous term in the Euler equation may then be eval- 
uated by use of spherical coordinates, 

r]Av + [^rj + C]'ViV-v) = ^f dr^-^drV^-^v . (15) 

Thus it is always the combination ^ = gjy + C that enters. 
Furthermore, it follows that a Hubble-type expansion, 
v{r) '^ r, is dissipation free in any dimension. 



B. Equations of motion 

The fluid-dynamic equations of motion reflect the con- 
servation of (baryon) charge, momentum, and energy. 
We are interested in the dynamics of small deviations 
from uniformity in a semi-infinite configuration and we 
focus on harmonic disturbances, 

p{r,t) =. po + Sp{x,t) = PO + Pfee^'"-'"' , (16) 
£(r,t) = eo + Seix,t) = eq + e-te''^""'"* , (17) 
p{r,t) = po + Sp{x,t) = po+PfcC^""-'"* , (18) 

and similarly for the other dynamical variables. 

The conservation of charge is ensured by the continuity 
equation, df^N'^ = 0, which here becomes 



C : dtp = -podxv =^ ujpk ^ pokvk 



(19) 



It serves to eliminate the flow velocity, Vk = iopk/{pok). 
The momentum equation simplifies considerably for the 
present scenario of small disturbances, 

M : hodtv = ~dx[p ~ Cdxv] ~ OxTTxx ~ dtq , (20) 

where ho = pq + Eq is the enthalpy density of the uniform 
system and the heat flow is q = (g, 0, 0) (see below). The 
equation for energy conservation is similarly simplified, 

E : dts ^ -hodxv - d^q ■ (21) 

By combining these latter two equations, (|20p and (HU, 
one obtains the sound equation. 



dtE-dxM: d^e = dJA[p ^ (dxv] + d^ 



(22) 



which amounts to oj'^Sk == k^pk — i£,{u!/ po)k'^pk where 

wp rei-a1l that f^ = in -i- t . see F,n (11 4-11. 



C. Dispersion equation 

When heat conductivity is ignored (k — 0), the energy 
density tracks the charge density, as follows immediately 
from the energy equation, p^ek = hopk- Furthermore, in 
the absence of a gradient term in the equation of state 
(see later), we have pk = Ps£k +PpPk with p^ = dePo{e, p) 
and Pp = dpPo{e,p) where po{£,p) is the microcanonical 
equation of state. Since the isentropic sound speed Vs 
is given by v'^ = Pe + (po/^o)Pp, we obtain the familiar 
viscous dispersion equation, w^ = v'^k'^ — i^{uj/ho)k'^ [3|. 

To obtain the dispersion equation with heat conductiv- 
ity included, we must invoke the form of the heat current, 

T up' 

q w -nldxT + Todtv] : qk = -iK[kTk —pk] ■ (23) 

Po k 

Insertion of this expression into the energy equation (j2ip 
yields a relationship between pk, £k, and Tk, 



£k 



ho k 


ho 


fc2 


— Pk + —qk - 


a —Pk ~ 


- IK Tk 


Po UJ 


PO 


UJ 



(24) 



where the term ~ npk has been ignored because if is 0{k) 
in comparison with {ho/ Po)Pk- Thermodynamics enables 
us to express 5e in terms of 5T and 5p, 



£k 



©/^^^(l)/^ 



-CvTk —pk , (25) 



where (Tje = d'^ao{e, p) and a^p = dedpao{s, p) are second 
derivatives of the entropy density (To(e, p). We have also 



used {de/dT),_ 



-l/Tase, the heat capacity at 



constant density. Using furthermore 



hoc^ee + Po<^p 



dp 
Ik 



(26) 



we may then obtain Tk from the energy equation (j2ip . 



Tk 



1 



To fdp 



I + iKk'^/cuCy Po \de 



Pk 



(27) 



The canonical equation of state pt{p) allows us to ex- 
press the pressure variation in terms of the variations in 
temperature and density, 

Pk=[l^]Tk+[Tr) Pk=\^] c,Tk + —VTPk, 



dT 



dp 



de 



Po 



(28) 
where vt is the isothermal sound speed, see Eq. ([5]). 

With these preparations, the dispersion equation can 
then be obtained by substituting the relations ([25]) , (|27p , 
(PSl) into the sound equation (E^ and using the relation- 
ship vl -vl = {T/h){dp/dT)p{dp/de)p, 



u 



2 7.2 



Vrpk 



ho 



9 2 



1 + ink^ /ujCv 



(29) 



retaining heat conduction terms only up to 0{k). This is 
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1. Gradient correction 

As noted above (end of Sect.|lT]), it is essential to take 
account of finite-range effects, without wliicli tlie spin- 
odal growtfi rate would become ever larger as the wave 
number is increased |f f| . Following Ref. [7[ , we introduce 
a gradient correction in the equation of state. To lead- 
ing order in the disturbance amplitudes, the effect of the 
gradient term on the local pressure is given by 



p{r) w _po(e(r-),p(r-)) -CpoVV(r) , 



(30) 



where po{e,p) is the microcanonical equation of state, 
i.e. the pressure in uniform matter having the specified 
energy and charge densities. The pressure amplitude is 
then modified accordingly, 

Pk ^ Pk + Cpok'^pk . (3f ) 

Hence we should augment the pk term in Eq. (I28|) . 

— VxPk ^ [ — Vt + Gpok^] Pk ■ (32) 

Po pa 

The full dispersion equation is then 



spinodal to the isothermal spinodal, vt{p,T) — 0, and 
generally increass the spinodal growth rates. 

For small distortions of uniform matter at given den- 
sity and temperature, the dispersion relation yields the 
eigenfrequencies Wk in terms of the equation of state 
Pt{p), the strength of the gradient term C, and the 
transport coefficients r]{p,T), C[p,T), and n{p,T). We 
describe below what we adopt for these key functions. 



IV. TRANSPORT COEFFICIENTS 

The deviation of the dynamical evolution from that of 
an ideal fluid is governed by three transport coefficients: 
the shear viscosity rj and the bulk viscosity C, (which here 
enter only through the effective viscosity ^ = |?7 -f C) as 
well as the heat conductivity k. Neither their magnitudes 
nor their dependencies on the environment (through p 
and T) are very well known. We shall therefore em- 
ploy simple parametrizations of their functional form and 
introduce one adjustable overall strength parameter for 
each one, thus enabling us to conveniently explore a range 
of physical scenarios. 



0.2 - 4fc2+c^fc^-z^ffc2+ ^ ^' ,;f fc^ (33) 

riQ riQ 1 + iK.k'^ /ijjCy 

i.e. the gradient term '^ Cfc'* is simply added, just as 
when there is no heat conductivity [7|. 



2. Solution of the dispersion equation 

When heat conductivity is included, k > 0, the disper- 
sion equation is of third order and, consequently, there 
are three eigenmodes for each wave vector fc, as one would 
generally expect since the energy density e is now no 
longer tied to the baryon density p. This equation has one 
purely imaginary solution, w^ = «7^, and a pair of gen- 
erally complex solutions, w^ , which are either both also 
imaginary, uJj^ = 17^. , or have the form w^ = =Fefc -I- 17^. 
In the latter case it is easy to see that 7fc < in the 
normal region where f^ > 0. 

The two frequencies wj obtained for ideal [i.e. non- 
dissipative) fluid dynamics are either purely real (out- 
side the isentropic spinodal phase region where v1 > 0) 
or purely imaginary (inside the isentropic spinodal re- 
gion where v'^ < 0). Thus, in ideal fluid dynamics, the 
region of spinodal instability is bounded by the isen- 
tropic spinodal, Vs{p,T) = 0. The introduction of viscos- 
ity adds a negative imaginary amount to the frequency. 



We then have 



Wfc 



I A vise ^2 to first order m 



^ = gT; -I- C, where we have introduced the character- 
istic viscous length Avisc(p, T) = £,{p,T)/ho{p,T)c. But 
the inclusion of viscosity does not change the region of 
instability. In contrast, the inclusion of heat conductiv- 

itv p.Ynandfi thp rpp'inn nf infit.j^.bilitv frnm thp ifipntrnnir 



A. Viscosity 

Using string theory methods, Kovtun et al. [12] made 
general arguments that the shear viscosity 77, in any rel- 
ativistic quantum field theory at finite temperature and 
zero chemical potential, has a lower bound, 77 > ?icr/47r, 
where a is the associated entropy density. 

It might therefore seem natural to use ry = rj^ha /At:, 
i.e. simply scale the minimum value by the factor rjQ >1. 
However, the entropy density a^ while appropriate in the 
context of ultrarelativistic nuclear collisions where the 
medium has a high temperature and a very small net 
baryon density, is not suitable in the present context 
where the focus is on matter at high baryon density and 
relatively modest temperature [13i] . A more appropriate 
quantity, suitable in both limits, is the enthalpy density 
/i(/3, T) = p + £ [ij]. When the net baron density p van- 
ishes it becomes h = Ta, whereas h/c^ approaches the 
mass density at high baryon density and low tempera- 
ture, h « mc^n ^ Ta, where n is density of particles 
and m is their mass. (For cold nuclear matter we have 
n K, p, since the pions and anti-nucleons have negligible 
populations, and hence m = rriN-) 

Furthermore, one would expect the viscosity to be pro- 
portional to the interparticle spacing d = 1/v}/^ [13J . 
which provides a convenient measure of the mean free 
path in a dense fiuid. When plasma has no net baryon 
density, d is inversely proportional to the temperature T, 
Tic/T = Ancod. The conversion constant cq is given by 



1 



Co 



(5s 



§59 



,C(3) 



0.12779 



(34) 



Therefore, based on these considerations, we shah make 
the following ansatz, 

V{P,T) ^ T]o—d{p,T)h{p,T) , r/o > 1 , (35) 

c 

which amounts to 77 = 77o?icr/47r in the baryon-free 
plasma, while it becomes 77 ~ r]QComcnd in the non- 
relativistic gas. This latter expression is similar to the 
familiar expression from classical transport theory for the 
shear viscosity coefficient in a dilute one-component gas, 
77 « ■^mvTi£, where v is the mean particle speed and £ is 
its mean free path. 

Since the bulk viscosity ( is generally expected to be 
significantly smaller than the shear viscosity 77, we shall 
take the effective viscosity to be ^ = I?/ + C ~ ^V- It 
is convenient to introduce the associated characteristic 
length which is proportional to the interparticle spacing. 



Kisc{p,T) 



1 ap,T) 

ch{p,T)/c^ 



VoCq d{p,T) . (36) 




Wave number k (fm" ) 



FIG. 2: The growth rate ykiPjT), as a function of the 
wave number k, calculated with finite-range fluid dynamics 
at p = 6.5 ps and T = 70MeV for four different combinations 
of dissipation: no dissipation (t^o = 0, kq = 0); minimal vis- 
cosity but no heat conduction (rjo —1,ko =0); no viscosity 
but minimal heat conduction (770 =0, kq =1); both minimal 
viscosity and minimal heat conduction {■qo = l,Ko —1). 



B. Heat conductivity 



The thermal conductivity is fundamentally related to 
the viscosity because it derives from the same microscopic 
transport processes. In a dilute classical gas it can be 
expressed as k « ^v£cv, where v is the mean particle 
speed and Cv = dT£T{p) is the heat capacity (equal to §71 
for a dilute classical gas). Since p = mv and h « mc^n, 
we see that K/77 w Cy/{h/c^). We therefore make the 
following ansatz. 



k{p,T) = KoCocd{p,T)cy{p,T) , kq > 1 



(37) 



with Co as given in Eq. p4p and assuming that the overall 
strength factor kq should be at least unity. Although the 
relation k/t] « c^/ {h/ c^) would suggest that the two nor- 
malization constants should be similar, kq ~ tjq, we pre- 
fer to leave them separately adjustable in order to make 
it possible to explore the effects of the two distinct types 
of dissipation. The characteristic length scale associated 
with heat conduction is then 



Kca.tiP,T) 



1 <P,T) 
cCv{p,T) 



Kq Co d{p, T) 



(38) 



While the approximate expressions for the transport 
coefficients, Eqs. (PSJ) and (1571) . should not be expected to 
be accurate, they will serve well for our present purpose 
of exploring the effect of the dissipative mechanisms on 
the spinodal decomposition, since they enable us to easily 
control the overall dissipative effects. 



Growth rates 



The spinodal isothermal and isentropic boundaries de- 
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to the thermodynamic limit of very long wave lengths, 
/c -T- 0. When the wave number k is increased, the region 
of instability steadily shrinks as the gradient term gives 
an ever larger contribution to the pressure. Thus, at a 
given temperature T, the lower spinodal boundary den- 
sity pA{k]T) increases steadily with fc, while the upper 
spinodal boundary density pB{k\T) decreases steadily. 
For a fixed value of A;, a contour plot of the growth rate 
7fc in the (p, T) phase plane therefore exhibits a ridge be- 
tween those two boundaries. The height of of the ridge 
decreases steadily as T is increased, first rather gently 
due to the dominance of the fermions at low tempera- 
tures. The local value of the maximum wave number for 
which instabilities exist, A:niax(/0, T), will have a similar 
appearance. 

We now illustrate, in Fig. [51 the resulting dispersion 
relations for thermodynamic scenarios relevant to the 
present study. Selecting a phase point in the central 
region of the phase coexistence region where both the 
isothermal and the isentropic sound velocities are imag- 
inary, p = Q.bps and T = 70MeV (see Fig. [1]), we con- 
sider the growth rate 7 as a function of the wave num- 
ber k of the density undulation being amplified. The 
non-dissipative treatment with ideal finite-range fluid dy- 
namics provides a convenient reference result. It yields a 
fastest growth time of about to ~ f .38fm/c which occurs 
for wave numbers near /cq ~ 2.0 fm~ , corresponding to 
an optimal wave length of £0 = 27r/fco ~ 3.14 fm. 

Relative to this reference, the inclusion of viscosity 
slows the growth but does not change the domain of in- 
stability which is still delineated by the vanishing of the 
isentropic sound speed Vs- We see that the inclusion of 
a minimal amount of viscosity (770 = 1) leads to a sig- 
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their strengths vary in approximate unison; the effect of 
varying their relative strengths may be judged from the 
results shown in Fig. [2j 

As expected from the discussion in the beginning of 
this section, the fastest growth at a given temperature oc- 
curs about midway between the corresponding spinodal 
boundaries. Furthermore, the temperature generally re- 
duces the growth rates, though only weakly at small T. 
The curves in the lower panel of Fig. [3] do not exhibit 
a monotonic decrease because the selected path along a 
constant density does not follow the ridgeline: the max- 
ima in the curves shown in the upper panel shift steadily 
towards smaller densities as T is increased, as would be 
expected from the general appearance of the phase dia- 
gram. 

Finally, we note that a five-fold increase in the dissipa- 
tion above the minimal value, i.e. changing (770, kq) from 
(1, 1) to (5,5), leads to a reduction in the growth rates 
7fc by only about a factor of two. 



FIG. 3: The maximum growth rate 70 (c/fm) shown at vari- 
ous compressions for fixed T — 70MeV {upper panel) and at 
various temperatures for fixed compression p/pa — 6.5 {lower 
panel), for three different assumptions about the dissipation, 
namely no dissipation: {rio,Ko) ~ (0,0); minimal dissipation: 
(1,1), and strong dissipation: (5,5). The open square cor- 
responds to the phase point at which the dispersion relation 
shown in Fig. [2] was obtained (using minimal dissipation). 



V. DYNAMICAL EVOLUTION 

After the above preparations, we are now in a position 
to address the dynamical evolution of the unstable collec- 
tive modes in the spinodal region of the phase diagram. 



scale towards larger values. 

On the other hand, relative to the ideal scenario, the 
inclusion of heat conductivity enlarges the domain of in- 
stability, the boundary being now determined by the van- 
ishing of the isothermal sound speed vt. Thus, generally, 
the inclusion of heat conductivity increases the growth 
rates, particularly at the high end of the unstable k range. 

While the inclusion of both minimal viscosity and min- 
imal heat conduction necessarily enlarges the unstable 
k range, it does somewhat reduce the fastest growth 
rates. However, it hardly affects the scale of the fastest- 
growing modes, /cmax- As the strengths of the dissipative 
terms are further increased, the growth rate 7^ decreases 
steadily and, at the same time, the maximum in the dis- 
persion relation moves gradually downwards in k. 

These features are present throughout the unstable re- 
gion of the {p, T) phase plane. In order to obtain an 
impression of how the growth rates change with p and T, 
we show in Fig. |3] the fastest growth rate 70 as we move 
away from the phase point explored in Fig. [2l either at 
constant T {upper panel) or constant p {lower panel), as 
indicated in Fig. [TJ In these illustrations we show the 
ideal scenario where no dissipation is present {rjQ, kq — 0), 
the scenario with minimal dissipation {rjQ, kq — 1), and a 
scenario with five times stronger dissipation {rjo , kq = 
thus covering the range suggested by the RHIC data 
llq : the currently favored estimate is rjo « 1-3. Since, as 
noted above, the heat conductivity is fundamentally re- 
lated tn fhp .shfiar visnnsit.v. it. .shnnid he PYnpct.fid that 



A. Dynamical phase trajectories 

We first specify dynamical phase trajectories, 
{p(t), T{t)), that are representative of the bulk matter in 
the collision zone, making use of the results presented in 
Rcf. 17]. In that work, a number of different dynamical 
models were used to extract the time evolution of the 
net baryon density, p{t), and the energy density, e{t), 
in the center of a head-on gold-gold collision for the 
range of collision energies anticipated at FAIR. The 
resulting dynamical trajectories in the (p, e) phase plane 
were remarkably independent of the specific model, in 
large part probably due to the robust nature of the 
mechanical densities which are subject to conservation 
laws, in contrast to the corresponding thermodynamical 
variables {ii,T). Nevertheless, there were significant 
variations in the detailed behavior. This is illustrated in 
Fig. m which shows density evolutions p{t) obtained with 
the 3-fluid model [ll| and UrQMD [li| at beam kinetic 
energies of 5 and 10 GeV/A. (Note that these beams 
are bombarded onto stationary targets, as will be done 
at FAIR; the same collision energies can be obtained in 
a collider configuration by using total beam energies of 
K, 1.8 and « 2.4 GeV/A, respectively.) We utilize these 
results to construct the dynamical phase trajectories 
considered. 

As expected, the degree of amplification achieved de- 
pends strongly on the length of time spent in the phase 
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FIG. 4: The time evolution of the net baryon density, p(t), 
at the center of a head-on gold-gold collision for bombarding 
energies of 5 and 10 GeV/A, as calculated with the 3-fluid 
[H] and the UrQMD O models (from Ref. 17]). 
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FIG. 5: Dynamical phase trajectories based on the 3-fluid 
and UrQMD density evolutions (shown in Fig.|4]) obtained for 
5 GeV/A in Ref. [13]; the associated time-dependent growth 
rates 7fc(t) are illustrated in Fig. [6] The symbols along the 
trajectories are equidistant in time with At — 1 fm/c, while 
the open dots on the left indicate the freezeout locations for 
bombarding energies of i? = 1, . . . , 10 GeV /A obtained from 
fits to experimental data as discussed in Ref. [201 ]. 



key feature, we consider in some detail the phase tra- 
jectories depicted on Fig. [5] The density evolutions p{t) 
are those calculated with the two models for 5 GeV/ A, 
while the temperature evolutions T{t) are (somewhat ar- 
bitrarily) prescribed with an eye towards the freeze-out 
conditions extracted from data ^20]. The essential dif- 
ference between these two trajectories is that UrQMD 
yields compressions that reach all the way into the decon- 
fined phase region, whereas the maximum compression 
achieved with the 3-fluid model lies inside the unstable 
region. Such a situation would be obtained with UrQMD 
as well at a somewhat lower collision energy (3-4 GeV /A) . 
When the turning point of the phase trajectory lies 



FIG. 6: The spinodal growth rate 7fe(t) — Ke[iJig{t)] for modes 
with k = 2fm~'^, calculated with minimal dissipation along 
the two dynamical (p, T) phase trajectories shown in Fig. [5| 



spinodal amplification for a longer time and, presumably, 
the resulting degree of amplification would be maximized. 
This key fact is illustrated in Fig. [51 For the penetrating 
(UrQMD) trajectory, the modes are exposed to amplifi- 
cation only during the two relatively brief periods when 
the phase point is traversing the spinodal region (and the 
amplification gained during the compressional traverse 
is largely lost during the high-density stage due to the 
equilibration process so only the amplification received 
during the expansion traverse is relevant). By contrast, 
the more optimal (3-fiuid) trajectory provides amplifica- 
tion for a sustained period of time, thus making a phase 
separation more likley to occur. 

While it is yet difficult to make specific predictions 
about the optimal collision energy, it seems evident that 
such an energy range exists, since the location of the 
(p, T) phase point associated with the maximum com- 
pression (the turning point) moves steadily downwards 
as the collision energy is lowered. The experimentally 
known freeze-out temperatures (indicated on the left in 
Fig. [5]) provide a lower bound on the phase region that 
could be accessed by nuclear collisions. The turning point 
is therefore likely to traverse the unstable phase region 
at fairly high temperatures where it is relatively nar- 
row. Consequently, the optimal range of collision energy 
is most likely not so wide. Because of this generic ex- 
pectation, we shall take the above phase trajectories as 
being reasonably representative of what may happen in 
the bulk region of an actual collision. 



B. Dissipative collective dynamics 

As the bulk of the evolving system traverses the un- 
stable phase region, its density fluctuations may be am- 
plified. In order to make quantitative estimates of this 
effect, we employ the method developed in Ref. [2lj for 
the evolution of collective modes subject to a dissipative 
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non-collective modes in the system. 

The treatment considers an ensemble of Af macroscop- 
ically similar systems {n} that have all been prepared 
with the same average density po relative to which indi- 
vidual perturbations are introduced, 



P^'^Kr) =. po + 6p(''>ir) 



in). 



Po 



fc#0 



(39) 



Generally, such as during the expansion stage of a nu- 
clear collision, the bulk density is time dependent. We 
shall assume that the effect of this overall evolution may 
be taken into account approximately by using the cor- 
responding time-dependent transport coefficients in the 
formulas below. We should then think of the mode in- 
dex fc as a mode number K rather than a wave number, 
since the expansion primarily stretches the modes with- 
out introducing much mode mixing [22|. In this way we 
may obtain the spatial size of the resulting fluctuations 
by scaling the initial wave length with the linear expan- 
sion factor, [po{ti)/po{tf)]^^'^, where D is the effective 
dimensionality of the expansion. 

The primary object of study is the associated spatial 
correlation function, cr(ri2) =~< Sp{ri) Sp{r2)* >-, where 

-< 5p{r,) Spir^r >^^Y. ^P^"^ (^i) '^/'^"^ (^2)* (40) 

n 

is the ensemble average and ri2 = ri — r2. The Fourier 
transform of cr(r) provides a convenient measure of the 
degree of density fluctuation at a given scale, since it is 
the ensemble average of |pfcp. 



= -< \Pk\ ^ = 



l^e-^V(r) 



(41) 



We are generally interested in the evolution of the col- 
lective modes in the system. Adopting the treatment of 
Ref. 12l[, we assume that the amplitude pk of a given 
collective mode k is governed by an equation of motion 
having the form 



-T;Pk{i) = -iujk{t)pk{t) + Bk{t) 
at 



(42) 



The complex eigenfrequencies uj^ = ek + ilk are deter- 
mined by the dispersion equation derived above, while 
the Brownian term Bk describes the residual coupling 
between the collective mode and the reservoir, i.e. the 
non-collective modes of the system. It is fluctuating in 
nature and is assumed to be Markovian, 



<Bk{t)Bk'{t'Y^^ 2Vkk'{t)S{t-t') 



(43) 



The equal-time correlation coefficient for two collective 
modes, akk'it) =^ Pk(t)pk'{t)* ^, is then given by [lH, 

^kk' it) - <Jkk' iU) e-'^kk' * + 2 J dt' Vkk' (f) e'^kk' (*'-*) , 



where cokk' = Wfc — w^/ . It is readily seen that it satisfies 
the following differential equation. 



dt 



O'kk'it) 



-iojkk' (i)Tfcfe' (t) + IVkk' (t) , (45) 



which was dubbed the Lalime equation in Ref. 2l|. 

For our present studies, we are particularly interested 
in the time evolution of the diagonal components of the 
covariance matrix, (Jkk — o'l^, which are equal to the fluc- 
tuation coefficients introduced in Eq. (|4T|) . Then 



'^lit) = 



al{t,)+ j\vk{t')e-^^ki^">dt' 



.2rfc(t) 



(46) 

where 'Dk = 'Dkk and we have introduced the following 
amplification coefficient, 



Tk{t) 



l-m[ujk{t'')\dt' 



jk{t')dt' . (47) 



If the environment is stationary, then 7^ and Vk re- 
main constant in time, so Tk{t) — 7fe(i — ti), and we 
obtain a simple exponential time evolution. 



<^lit) 



Vk_ 

7fc 



,27fc(t-t,) 



- 1 



= (iOe^ 



27fc(t-*i) 



(48) 



When 7fc is negative, as is normally the case, the mode is 
stable and relaxes exponentially towards its equilibrium 
value al, <Tk{t) — > —Vk/jk, the associated relaxation 
time being tk = l/7fc- More generally, within the stable 
regime, <T^(t) will seek to relax towards its instantaneous 
equilibrium value fj|(i) = — X'fc(i)/7fc(i), in accordance 
with the Einstein relation. In the opposite case, when 7^ 
is positive, the mode is unstable and the fluctuations ex- 
hibit an exponential growth, with both the original mode 
fluctuations (Jkitt) and the noise Vk/lk being amplified. 
The diffusion coefficients Vk represent the coupling 
of the collective modes to the residual system. They 
thus play two important roles in the dynamical evolu- 
tion of the density fluctuations. In the stable regime (as 
mentioned above) the coupling determines the relaxation 
times tk for the relaxation of the fluctuation coefficients 
(7^ towards their appropriate equilibrium values (t| , while 
in the unstable regime it continually produces additional 
fluctuations that are also amplified. Thus, even in the 
absense of initial fiuctuations, the coupling of the col- 
lective mode to the residual system will continually cre- 
ate fluctuations that will subsequently become ampliflcd 
or damped, as governed by the Lalime equation (j45l) . 
In particular, the diffusion coefficients enable the sta- 
ble modes to continually adjust their fluctuations as the 
equilibrium variances evolve. 



C. Onset of phase separation 



The amplification coefficient Tk defined in Eq. (|47l) is 
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cation factor, 



Gk = exp( / -fk{t)dt) 



(49) 



Wave number k (fm" 



FIG. 7: As a function of the wave number k is shown the 
ampUfication factor Gk (see Eq. (|49|l 'l resuhing from motion 
along the 3-fluid phase trajectory displayed in Fig. [S] for var- 
ious degrees of dissipation (indicated by the values of ryo and 
Ko). Also shown is the result for the UrQMD trajectory in 
Fig. [5]using minimal dissipation, i.e. {rjo, kq) — (1, 1). 



interval considered. When there is no dissipation, cj^ is 
always real so jk vanishes outside the unstable region. 
But in the presence of dissipation, jk is negative out- 
side the unstable region, and increasingly so the stronger 
the dissipation. As a result, the local relaxation time 
is reduced which in turn ensures that the fluctuations 
remain close to their equilibrium value until shortly be- 
fore the trajectory enters the unstable region. However, 
by the same token, any amplification acquired while the 
trajectory is inside the unstable region is correspondingly 
quickly lost after the trajectory has reentered the normal 
region and again equilibrates. In fact, the dissipative 
length tends to be larger at late times, due to the in- 
creased particle spacing in the more dilute system, thus 
shortening the relaxation time. 

Our present study aims to clarify the prospects for 
the spinodal amplification to become sufficiently large to 
cause a phase-separating clumping of the system. Since 
our treatment is perturbative, it is only valid as long as 
the fluctuations remain relatively small. The occurrence 
of large fluctuations in the calculation should then be 
taken as a signal that clumping is likely. It should be 
recalled that the bulk of the system is thermodynam- 
cially unstable so even a modest degree of fluctuation 
may cause a catastrophic breakup. The prospects for 
this to occur are enhanced by the fact that fluctuations 
created in the mechanically unstable (spinodal) region 
may become futher amplified during the traverse of the 
adjacent metastable region, just as impurities introduced 
in this region may trigger condensation. 

Thus the key question is how much amplification may 
occur during the unstable era. To elucidate this issue in 
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where the integral is only over those times during which 
the mode is unstable, 7^ > 0. 

Figure [7] shows values of G^ obtained for the entire 
range of wave numbers, as obtained with various degrees 
of dissipation. While we concentrate on the phase tra- 
jectory that reaches its maximum compression inside the 
spinodal region (the one based on the 3-fluid model re- 
sults for 5 GeV/A) , we also show the result of a more pen- 
etrating trajctory to bring out the importance of tuning 
the collision energy for optimal effect. 

In addition to the result of ideal fluid dynamics, we 
show results for various degrees of dissipation, ranging 
from minimal, i.e. (770,^0) = (lil)j to five times that. 
While viscosity generally slows the evolution, thus also 
supressing the growth of instabilities, heat conductivity 
generally increases the growth rate. The combined effect 
of introducing small amounts of dissipation then tends to 
enhance the amplification. Thus the resulting degree of 
non-uniformity is fairly robust against moderate changes 
in the dissipation strength and, consequently, our con- 
clusions do not appear to be sensitive to the specific 
parametrizations of the transport coefficients. 

The results displayed in Fig. [7] bring out the charac- 
teristic feature of spinodal instability, namely that the 
amplification mechanism favors certain length scales. We 
note that the two-point correlation coefficient af. is pro- 
portional to G\ and thus exhibits a stronger peaking. 
More generally, since the iV-point correlation is propor- 
tional to G^ , the spinodal effect manifests itself progres- 
sively stronger in the higher-order correlations. 

We see that Gk is peaked around /cq w 1.6-1.3 fm~^, 
depending on the degree of dissipation, which corre- 
sponds to wave lengths Iq — 27r/fco ~ 3.9-4.8 fm. If, at a 
temperature of 80-100 MeV, a uniform system situated 
at the lower edge of the spinodal region (hence having 
a density po ~ 5ps) transforms itself into plasma drops 
embedded in a hadron gas, with each subsystem having 
the corresponding coexistence density, pi{T) w ips and 
P2{T) ~ 8ps, then we may expect the drop radius to be 
i?drop « ^Uipo - Pi)/(P2 - Po)]'/' « 1.2-1.5 fm. Each 
such drop would then have a baryon number of 10-18. 
Since this is only a relatively small fraction of the col- 
lision system one would expect that several such drops 
would be formed during the phase separation. On the 
other hand, a drop of such a size is sufficiently large to 
be regarded as a macroscopic source that will undergo 
statistical hadronization. 



VI. CONCLUDING REMARKS 

In the exploration of the phase diagram of strongly in- 
teracting matter by means of nuclear collisions, the mech- 
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unique signals of the first-order phase transition. For the 
planning of possible experimental campaigns to search 
for spinodal phase separation, it is useful to have esti- 
mates of which bombarding energies are expected be be 
optimal for producing the effect. The analysis presented 
above suggests that these bombarding energies lie at the 
lower end of the anticipated FAIR range. It is also well 
within the energy region proposed for NICA but appears 
to be too low for experiments at RHIC to be feasible. 

However, when making a quantitative prediction of the 
optimal collision energy range, it should be kept in mind 
that the calculated dynamical phase trajectories [17], on 
the basis of which we made our estimates, extracted the 
conditions right at the center of the collision zone. There- 
fore they most likely provide an upper limit on the com- 
pression and it must be exptected that the average com- 
pression and excitation over an extended volume will be 
smaller than those shown in Figs. |3] and [S] Consequently, 
the optimal bombarding energy is probably somewhat 
higher than what the above idealized analysis would sug- 
gest. For example, it may well be that the 5 GeV/A 
UrQMD simulation, whose maximum compression over- 
shoots the spinodal region, would be quite suitable for 
spinodal clumping because the average degree of com- 
pression over an extended region favors phase separation. 

It should also be recognized that the spinodal phase re- 
gion is surrounded by a thermodynamically metastable 
region through which the phase trajectory has to pass 
after the spinodal mechanism has acted. During this tra- 
verse, sufBciently significant deviations from uniformity 
will be further amplified and phase separation is thus 
more likely to occur than suggested by just the pertur- 
bative treatment employed here. In fact, the spinodal 
amplification of the preexisting fiuctuations may be re- 
garded as merely providing seeds for a subsequent non- 
linear breakup evolution in the metastable region. 

We note the analogy with the search for spinodal frag- 
mentation [5] as a signal of the nuclear liquid-gas phase 
transition jS]. The collision energy had to be carefully 
adjusted to producing the phenomenon: a too high en- 
ergy would yield an explosive vaporization, while a too 
low energy would cause the two nuclei to fuse (and sub- 
sequently deexcite by light-particle emission). But in the 
optimal energy range, the bulk of the combined system 
would first compress (to 2-3 times normal) and then start 
to expand; however, the expansion would tend to stall, 
leaving the system in a dilute (and nearly spherical) con- 
figuration for a sufficient length of time (several times the 
typical growth time) to cause the most unstable modes 
to grow dominant, thereby leading the system towards a 
breakup into nearly equal-sized intermediate-mass frag- 
ments. Our present analysis suggests that collisions at 
suitably tuned relativistic energies would also cause the 
bulk of the system to spend several growth times inside 
the unstable phase region and thus enable the spinodal 
formation of plasma drops. 

In the liquid-gas case the production of equal-size nu- 
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ambigous signal of the spinodal breakup mechanism and 
thus for the existence of a first-order transition. The 
confinement transition is inherently more difficult to in- 
vestigate experimentally because any plasma drops that 
may have been formed will ultimately hadronize and are 
thus harder to identify. Nevertheless, the transient ex- 
istence of such spatially separated blobs of deconfined 
matter may be revealed by careful examination of suit- 
able multi-particle correlations. While some relatively 
schematic studies have already been made for the pur- 
pose of identifying such observables J23l - l25| there is a 
need for much more refined studies. We hope that the 
present investigation, which suggests that spinodal phase 
separation might indeed occur, will provide an incentive 
for such endeavors. 
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Appendix A: Equation of state 

The present study requires an equation of state of 
strongly interacting matter that displays the expected 
phase structure. Although significant progress has been 
made in understanding the thermodynamical properties 
of each of the phases separately, our current understand- 
ing of the phase coexistence region is not yet on firm 
ground. We therefore employ a conceptually simple ap- 
proximate equation of state which will suffice for our 
present explorations. 

For this purpose, we approximate the confined phase 
by an ideal gas of nucleons and pions augmented by a 
density-dependent interaction energy, while the decon- 
fined phase is taken as an ideal gas of quarks and gluons 
with their interactions described by a bag constant. The 
desired phase structure is then generated by suitable in- 
terpolation between these two pure phases. 



1. Confined phase 

The confined phase is approximated as an ideal gas 
of pions, nucleons, and antinucleons, plus an interaction 
term. The total hadronic pressure is thus 

p" =Pit+Pn +Pn +Pw, (Al) 

where the contribtion from the ideal pion gas is 

P7r[T) = -g^T / 
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FIG. 8: The equation of state, pt{p), for the (p,T) region 
relevant for ordinary nuclear matter: the dependence of the 
pressure on the compression for specified T (shown in MeV). 



with (^TT = 3 and 771,^ — 140 MeV, while the nucleons and 
antinucleons contribute 
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respectively, with gjv = 2 x 2 = 4 and rriM — 940 MeV. 
The net baryon density p^ = p^ — ppf then follows, 
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pede sinh /3/io 
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(A5) 



and the entropy density is a" — dp" / dT . Finally, the 
contribution from the interaction energy density 'w{p) is 
Pw{p) = pdpw(p) — w{p) and the parameter po is related 
to the chemical potential p hy p = ptQ + dpW. 

The interaction energy density w(p) has the form 



w{p) 
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where we use a — 1 and j3 = 2. The strength coeffi- 
cients A and B are then adjusted so that nuclear mat- 
ter saturates at ps = 0.153 fm~ and the associated com- 
pression modulus is, K = Km + K^, = 300 MeV, where 
Kn — —^Ep sa —43 MeV is the contribution from the 
Fermi motion of the nucleons (which is negative). The 
binding energy of nuclear matter is then also roughly re- 
produced. The resulting equation of state for nuclear 
matter is shown in Fig. [51 



2. Deconfined phase 

The deconfined phase is taken as an ideal gas of mass- 
less gluons and light quarks with a standard bag constant. 
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FIG. 9: Phase crossing: The pressures in the two ideal- 
ized phases are shown as functions of the chemical poten- 
tial /i for T = 0; the systems are in mutual thermody- 
namic equilibrium at the p value for which the two curves 
cross. The crossing points obtained by the same procedure 
for r > are connected by the solid curve which terminates 
at Tmax ~ 156.4 MeV. 



where pg = (7g(7r^/90)r^ with gg = 2x8 = 16 is the gluon 
pressure while the quarks and antiquarks contribute 
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with g^ = 2 X 3 X 2 = 12 and /ig = jp. The net baryon 
density in the plasma is then 



P^ = 



9pQ 



= nAiT" 



dp. 9 

while the entropy density is 
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For the bag constant we use B — 300MeV/fm'^. 



3. Interpolation 

At zero temperature and zero chemical potential, the 
pressure of the nucleon gas vanishes while that of the 
quark gas is equal to —B. The confined phase is then the 
thermodynamically favored one. However, as the chem- 
ical potential is raised, the plasma pressure increases 
faster than the hadronic pressure, so the two curves, 
p" (T — 0,p) and p'^{T — 0,p) cross at a certain value 
of p, above which the deconfined phase is favored, as il- 
lustrated in Fig. [3 This phase crossing procedure can be 
repeated for any temperature up to T^ax and the result- 
ing crossing points are included in Fig. [HI 

For the discussion of phase coexistence, it is conve- 
nient to work in the canonical representation where the 
temperature is specified. Then the condition of phase co- 

p."X"ifitpncp 7.. p.. sa.mp tpmnpratnrp chpmicf^l "nntpntij^.l and 
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FIG. 10: Phase crossing: For the range of temperatures 
where the two ideaUzed phases coexists, < T < Tmax ~ 
156.4 MeV, the associated coexistence densities are shown 
(dashed curves). The soUd curves show the corresponding 
coexistence densities for the spline-based equation of state. 



pressure at two different densities pi and p2, amounts 
to the condition that frip), the free energy density as 
a function of density, have common tangents. This is 
readily seen since prip) = dp frip) implies that the two 
chemical potential are then equal, /^i = /i2, and since the 
tangent at p^ is given by ti[p) = f^{p) + {p- Pi)fi{p) the 
fact that ti(p) = i2{p) immediately implies pi = P2- 

Figure [TT] shows fx^oip) and frf=oip)- The former 
curve starts at zero but grows more rapidly than the 
latter which starts at J5, so the two curves cross and, 
furthermore, since they both have positive curvature, a 
common tangent exists. Because of this generic feature, 
it is expected that cold matter exhibits a first-order phase 
transition when compressed. 

While the simple gas models employed presumably 
provide reasonable (though still somewhat idealized) de- 
scriptions of the two individual phases well away from the 
coexistence region, neither one is suitable in the phase- 
coexistence region. In order to describe the transition 
region, we represent the free energy density there by a 
S^'^-order polynomial /t(/o) that matches the values of 
It 7 dpfx 1 and 9p/|^ at a density p^ < p^ and the val- 
ues of /^, dpf^ , and dpf^ at a density p^ It Pt^ where 
/5y and p^ are those densities at which the two ideal- 
ized curves fxip) and frp{p) have a common tangent. 
For lower densities, p < pif , we use the idealized hadron 
gas, f^{p), and at higher energies, p > p^, we use the 
idealized plasma, frf{p). 

When the lower matching density p^ is chosen suffi- 
ciently close to pif and the higher matching density p^ 
is chosen sufficiently close to p^ then the resulting spline 
function frip) also has a common tangent and so the 
system has a first-order transition at the particular tem- 
perature T. The associated densities (where the common 

ta.no'PTit tnnrhps frT^( n\\ arp thpn tlip cnPYistpncp Hpnsitips. 
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FIG. 11: The free energy density at zero temperature, 
/y^o(p), as a function of the degree of compression p/pc- 
The two individual phases fo {p) and fffip) are shown to- 
gether with the connection between them, /o(p), obtained by 
sphning between the two open squares located at p^ and /5q . 
The resulting coexistence points (filled circles) are located at 
pi and P2; they are connected by the associated common tan- 
gent (the Maxwell line). The spinodal boundary densities pA 
and pB are indicated by the open circles. 



pi and p2, at that temperature. 

This spline procedure is illustrated in Fig. [TT] for T — 
and the corresponding pressure is shown in Fig. [T^l By 
suitable adjustment of the matching densities p^ and p^, 
it is thus possible to design an equation of state having 
the desired phase structure, namely a first-order phase 
transition that becomes ever weaker as the temperature 
is raised and terminates in a critical point at a finite 
density pc- The resulting phase coexistence boundaries, 
(pi,r) and (p2,T), are shown in Fig. [TUl the full phase 
diagram having already been shown in Fig. [T] 



Appendix B: Spline method 

We describe here a convenient spline method that en- 
ables us to match values and derivatives at two points. 

We seek an polynomial expression for the function f{x) 
(here the free energy density frip)) that matches the 
specified values ag = f{x ^ a) and bo = f{x — b) as 
well as the associated derivatives up to any order n, Oj = 



f(») 



{d'f/dx% 



and bj 



f, 



(») 



{d'-f/dx'')x=b, i 



Ja 

l,...,n. 

If fn-i{x) denotes the spline approximation that 
matches derivatives up to the (n — 1)"^ order (which is 
a polynomial of order 2?i — 1), then we may obtain the 
approximation for the next order n by writing 

fn[x) = fn-lix) + {b-x) [X-a] r . 

b — a 

(Bl) 
Since the factor {b — x)"(x — a)" and its derivatives up 

tn nrdpT n. — 1 va.ni.'sli at n and h it fnllnwrs that f„ (cr) 
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FIG. 12: The pressure at zero temperature, pT=o(p), as a 
function of the degree of compression p/pc, corresponding to 
the free-energy density shown in Fig. 1111 The two individual 
phases are shown together with the connection between them 
obtained by sphning between the two open squares. The re- 
sulting coexistence points (open circles) are connected by the 
associated Maxwell line. The spinodal boundaries are situ- 
ated at the two extrema (open circles). 



satisfies the matching conditions up to order n — 1, i.e. 
.fi'\a) = f^M = a, and f^\b) = f^,(^ = 6. for 
i < n. 

The only remaining task is thus to determine the two 
coefficients q;„ and /3„ which can be done by matching 



also the n derivative at a and 6, namely fa 



(n) 



(n) 

fh ^ bn- It is elementary to show that they are given 
by the following expressions, 



an~/l-i(a) 
n\{b--aY 



, Pn 



flC\{h) 



n\{a-bY 



(B2) 



It is possible to also determine the derivatives of the 
spline function by iteration. In particular, the derivatives 
at the two matching points are given by 






(B3) 



{n - v)\ {2v ~ n + 1)\ 
&Hb) = f!)-{ib) + 

{n - i^y. {2iy - n + ly. 



[{ly + l)a^ - (n - p)/3^] 



(B4) 



[{iy + l)/3^-{n-iy)a^] 



and 



where the first terms are present only for n < 2t/, while 
the second terms are present only for ly < n < 2^ + 1. 
We note that 

/(2"+i)(a) = (2n + l)!(-)"[/3„-a„] = fi'^+^Hb) , (B5) 

consistent with the fact that the highest non-zero deriva- 
tive is a constant. 

The above expressions can be used iteratively to de- 
termine the spline polynomial fn{x) for any oder n > 
as well as all of its 2n + 1 derivatives. 
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